///////////////////////////////////////////////////////////////////////////////
//
// labFisOndeSup.cc - acquire data for SNS Lab Fis Water Surface Waves exp.
// 
// 1st release: A. Lusiani <a.lusiani@sns.it> March 2004
//
// this program includes parts from poll.c, Part of Comedilib,
// Copyright (c) 1999,2000 David A. Schleef <ds@schleef.org>
// poll.c may be freely modified, distributed, and combined with
// other software, as long as proper attribution is given in the
// source code.
//

///////////////////////////////////////////////////////////////////////////////
//
// Modifications:
//
//

///////////////////////////////////////////////////////////////////////////////
//
// headers
//

extern "C" {
#include <stdio.h>
#include <stdlib.h>
#include <termios.h>
#include <comedilib.h>
#include <fcntl.h>
#include <unistd.h>
#include <errno.h>
#include <ctype.h>
#include <string.h>
#include <math.h>
#include <fcntl.h>
#include <sys/select.h>
}

#ifndef __CINT__
#include "TFile.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TRandom.h"
#include "TMath.h"
#endif

#include "labFisDaq.h"

///////////////////////////////////////////////////////////////////////////////
//
// function declarations
//

///////////////////////////////////////////////////////////////////////////////
//
// variable  declarations
//

const char* 		COMEDI_DEV = 		"/dev/comedi0";
int			AI_SUBDEV;
int			AO_SUBDEV;

//--- buffer size for streaming data out of the driver
const unsigned int	COMEDI_BUFSIZ =		(4*1024*1024);

//
// input channel & scan definitions
// - scan a total of AISCAN_NUM AI channels beginning with AISCAN_BEG
//

//--- Analog Input definitions for the two photo detectors
const unsigned int 	PHDET1_AI_CHAN =	0;
const unsigned int 	PHDET2_AI_CHAN =	1;

//--- Analog Input voltage range
const float 		AI_VMIN = 		-1;
const float 		AI_VMAX =		1;

//--- select analog in voltage reference from GROUND
const unsigned int	AI_AREF =		AREF_GROUND;

const unsigned int	AISCAN_BEG =		0;
#define			AISCAN_NUM		(2)
unsigned int		aiScanList[AISCAN_NUM];

//--- scan rate and buffer size
const unsigned int 	AISCAN_RATE = 		5000;
#define 		AISCAN_SAMPLES 		(30000)
sampl_t			aiScanData[AISCAN_SAMPLES*AISCAN_NUM];

//
// output channel & scan definitions
// - scan a total of AISCAN_NUM AI channels beginning with AISCAN_BEG
//

//--- Analog Out definition
const unsigned int 	TWEETER_AO_CHAN = 	0;

const double		TWEETER_AMPLITUDE = 	5;
//--- min voltage for driving tweeter
const double 		TWEETER_MIN =		-TWEETER_AMPLITUDE;
//--- max voltage for driving tweeter
const double		TWEETER_MAX =		TWEETER_AMPLITUDE;
//--- zero voltage for driving tweeter
const double		TWEETER_ZERO =		0;

//--- select Analog Output voltage reference from GROUND
const unsigned int	AO_AREF = 		AREF_GROUND;

#define			AOSCAN_NUM		(1)
unsigned int		aoScanList[AOSCAN_NUM];

//
// Analog output waveform generation
// (tested on Measurement Computing card)
//
// Limitations: 
// - analog output rate is limited to 10kHz (0.1ms)
// - one must preload at least 5000 samples
//
// output one full period of a sine wave with T = 0.1s
// - this corresponds to 1000 samples
// - must complete the buffer with 4000 samples of silence
// 

//--- duration of AO wave in ns
#define			AOWAVE_DURATION	(50000000)
#define 		AOWAVE_SAMPLES 	(1000)
#define 		AOSCAN_SAMPLES 	(5000)
const unsigned int	AOSCAN_TIME = 	AOWAVE_DURATION/AOWAVE_SAMPLES;
sampl_t 		aoScanData[AOSCAN_NUM*AOSCAN_SAMPLES];

const double		AOSCAN_FREQ = 1;

///////////////////////////////////////////////////////////////////////////////
//
// code
//

//
// comedi_write_wave
//
// write comedi device with the content of a waveform buffer
//
int comedi_write_wave(comedi_t *cifh,
		      sampl_t *waveData,
		      unsigned int waveDataNum,
		      unsigned int& waveDataCycle,
		      unsigned int numSamples)
{
  char* waveByte = (char*)waveData;
  size_t waveByteNum = waveDataNum*sizeof(sampl_t);
  size_t waveByteCycle = waveDataCycle*sizeof(sampl_t);
  size_t numBytes = numSamples*sizeof(sampl_t);

  size_t writeBytes;
  size_t writtenBytes;
  size_t availBytes;

  for(;numBytes>0;) {
    //--- total n. of bytes to write
    writeBytes = numBytes;
    //--- available bytes from aoScanCycle to end of waveform buffer
    availBytes = waveByteNum - waveByteCycle;
    if (availBytes<writeBytes) {
      writeBytes = availBytes;
    }
    writtenBytes = write(comedi_fileno(cifh),
			 &waveByte[waveByteCycle],
			 writeBytes);
    if (writtenBytes<0){
      return writtenBytes;
    }
    printf("ao scan written = %d %d\n", writtenBytes, waveByteCycle);
    numBytes -= writtenBytes;
    waveByteCycle += writtenBytes;
    if (waveByteCycle>=waveByteNum) {
      waveByteCycle -= waveByteNum;
    }
  }
  waveDataCycle = waveByteCycle/sizeof(sampl_t);
  return 0;
}

int main(int argc, char *argv[])
{
  int i;
  int rc;
  comedi_t *cifh;
  comedi_cmd cmd;
  comedi_cmd cmdAo;
  lsampl_t maxDataAi;
  lsampl_t maxDataAo;
  double minVoltsAi, maxVoltsAi;
  double minVoltsAo, maxVoltsAo;
  unsigned int rangetAi;
  comedi_range* rangepAi;
  unsigned int rangetAo;
  comedi_range* rangepAo;
  lsampl_t tweeterDataMin;
  lsampl_t tweeterDataMax;
  lsampl_t tweeterDataZero;
  //--- buffer for ni_pcimio Analog Input driver (16-bit unsigned samples)
  int pendingNewline = 0;
  unsigned int len;
  unsigned int scansTot;
  
  //--- endpoints of acquisition are valid numbers
  rc = comedi_set_global_oor_behavior(COMEDI_OOR_NUMBER);

  cifh = comedi_open(COMEDI_DEV);
  if(!cifh){
    perror(COMEDI_DEV);
    exit(1);
  }

  //--- set comedi device streaming buffer size
  rc = comedi_set_buffer_size(cifh, AI_SUBDEV, COMEDI_BUFSIZ);
  if (rc<0) {
    fprintf(stderr, "comedi_set_buffer_size: %s\n",strerror(errno));
    exit(1);
  } else {
    rc = comedi_get_buffer_size(cifh, AI_SUBDEV);
    if (rc<0) {
      perror("comedi_get_buffer_size");
    } else {
      fprintf(stderr, "comedi buffer size for AI set to 0x%x bytes\n", rc);
    }
  }

  //--- find which subdevice offers Analog Input
  rc = comedi_find_subdevice_by_type(cifh, COMEDI_SUBD_AI, 0);
  if (rc<0){
    comedi_perror("comedi_find_device AI");
    exit(1);
  }
  AI_SUBDEV = rc;

  //--- find which subdevice offers Analog Output
  rc = comedi_find_subdevice_by_type(cifh, COMEDI_SUBD_AO, 0);
  if (rc<0){
    comedi_perror("comedi_find_device AO");
    exit(1);
  }
  AO_SUBDEV = rc;

  //--- find AI range for -10 Volt .. 10 Volt
  rc = comedi_find_range(cifh, AI_SUBDEV,
			 PHDET1_AI_CHAN, UNIT_volt, AI_VMIN, AI_VMAX);
  if (rc<0){
    comedi_perror("comedi_find_range");
    exit(1);
  }
  rangetAi = rc;

  //--- get ptr to range struct
  rangepAi = comedi_get_range(cifh, AI_SUBDEV, PHDET1_AI_CHAN, rangetAi);
  //--- get maximum possible sample for AI
  maxDataAi = comedi_get_maxdata(cifh, AI_SUBDEV, PHDET1_AI_CHAN);
  minVoltsAi = comedi_to_phys(0, rangepAi, maxDataAi);
  maxVoltsAi = comedi_to_phys(maxDataAi, rangepAi, maxDataAi);

  fprintf(stderr, "Comedi offers AI range n. %d, %f .. %f\n",
	 rangetAi, minVoltsAi, maxVoltsAi);

  //--- build channel list for DAQ scan of two photo detectors
  for(i=0;i<AISCAN_NUM;i++){
    aiScanList[i] = CR_PACK(AISCAN_BEG+i, rangetAi, AI_AREF);
  }

  //--- find AO range for -10 Volt .. 10 Volt
  rc = comedi_find_range(cifh, AO_SUBDEV,
			 TWEETER_AO_CHAN, UNIT_volt,
			 TWEETER_MIN, TWEETER_MAX);
  if (rc<0){
    comedi_perror("comedi_find_range");
    exit(1);
  }
  rangetAo = rc;

  //--- get ptr to range struct
  rangepAo = comedi_get_range(cifh, AO_SUBDEV, TWEETER_AO_CHAN, rangetAo);
  //--- get maximum possible sample for AO
  maxDataAo = comedi_get_maxdata(cifh, AO_SUBDEV, TWEETER_AO_CHAN);
  minVoltsAo = comedi_to_phys(0, rangepAo, maxDataAo);
  maxVoltsAo = comedi_to_phys(maxDataAo, rangepAo, maxDataAo);

  fprintf(stderr, "Comedi offers AO range n. %d, %f .. %f\n",
	 rangetAo, minVoltsAo, maxVoltsAo);

  tweeterDataMin = comedi_from_phys(TWEETER_MIN, rangepAo, maxDataAo);
  tweeterDataMax = comedi_from_phys(TWEETER_MAX, rangepAo, maxDataAo);
  tweeterDataZero = comedi_from_phys(TWEETER_ZERO, rangepAo, maxDataAo);

  //
  // now prepare command for analog output
  //
  memset(&cmdAo, 0, sizeof(cmdAo));
  cmdAo.subdev = AO_SUBDEV;
  cmdAo.flags = 0;
  cmdAo.start_src = TRIG_INT;
  cmdAo.start_arg = 0;
  cmdAo.scan_begin_src = TRIG_TIMER;
  cmdAo.scan_begin_arg = AOSCAN_TIME;
  cmdAo.convert_src = TRIG_NOW;
  cmdAo.convert_arg = 0;
  cmdAo.scan_end_src = TRIG_COUNT;
  cmdAo.scan_end_arg = AOSCAN_NUM;
  cmdAo.stop_src = TRIG_COUNT;
  cmdAo.stop_arg = AOSCAN_SAMPLES;

  cmdAo.stop_src = TRIG_NONE;
  cmdAo.stop_arg = 0;

  cmdAo.chanlist = aoScanList;
  cmdAo.chanlist_len = AOSCAN_NUM;
  
  //--- assemble channel description for 1 analog output channel
  aoScanList[0] = CR_PACK(TWEETER_AO_CHAN, rangetAo, AO_AREF);

#undef CONTINUOUS_TEST
#ifdef CONTINUOUS_TEST
  //--- waveform: continuous sine wave
  for( i=0; i<AOSCAN_SAMPLES; i++) {
    //--- x varies from 0 to AOSCAN_SAMPLES/AOWAVE_SAMPLES
    double x = double(i)/(AOWAVE_SAMPLES);
    aoScanData[i] = comedi_from_phys(TWEETER_AMPLITUDE*sin(AOSCAN_FREQ*2*M_PI*x),
				     rangepAo, maxDataAo);
    // printf("x=%f, val=%d, i=%d\n", x, aoScanData[i], i);
  }
#else
  //--- waveform: one sine semi-period followed by silence
  for( i=0; i<AOWAVE_SAMPLES; i++) {
    //--- x varies from 0 to 1
    double x = double(i)/(AOWAVE_SAMPLES);
    aoScanData[i] = comedi_from_phys(TWEETER_AMPLITUDE*
				     // -pow(sin(AOSCAN_FREQ*1*M_PI*x),4)
				     -sin(AOSCAN_FREQ*2*M_PI*x)
				     ,
				     rangepAo, maxDataAo);
    // printf("x=%f, val=%d, i=%d\n", x, aoScanData[i], i);
  }
  for( ; i<AOSCAN_SAMPLES; i++) {
    aoScanData[i] = comedi_from_phys(TWEETER_ZERO, rangepAo, maxDataAo);
  }
#endif
  
  //--- check if output scan command is possible
  if (labFisScanCmdTest(cifh, &cmdAo)) {
    exit(1);
  }

  //--- set tweeter at TWEETER_ZERO
  rc = comedi_data_write(cifh, AO_SUBDEV, TWEETER_AO_CHAN,
			 rangetAo, AO_AREF, tweeterDataZero);
  if (rc<0){
    comedi_perror("comedi_write_data 1");
    exit(1);
  }
  //--- wait 1s, let water settle in case tweeter was not already there
  usleep(1000000);

  //--- prepare AI scan for photo-detector input
  labFisScanPrepare(cifh, AI_SUBDEV, &cmd,
		    TRIG_NOW, AISCAN_RATE, TRIG_COUNT, AISCAN_SAMPLES,
		    aiScanList, AISCAN_NUM);
  
  //--- start comedi Analog Input scan
  rc = comedi_command(cifh, &cmd);
  fprintf(stderr,"Start DAQ, rc= %d\n", rc);
  if (rc<0) {
    fprintf(stderr,"errno=%d\n",errno);
    comedi_perror("comedi_command");
    exit(1);
  }

  //--- wait 0.2s, let water settle in case tweeter was not already there
  usleep(200000);

#undef SQUARE_STIMULATION
#ifdef SQUARE_STIMULATION

  //
  // stimulate tweeter with square signal
  //

  for(int iii=0; iii<1; iii++) {
  //--- set tweeter at TWEETER_MAX
  rc = comedi_data_write(cifh, AO_SUBDEV, TWEETER_AO_CHAN,
			 rangetAo, AO_AREF, tweeterDataMax);
  if (rc<0){
    comedi_perror("comedi_write_data 2");
    exit(1);
  }
  //--- wait 25ms
  usleep(25000);

  //--- set tweeter at TWEETER_MIN
  rc = comedi_data_write(cifh, AO_SUBDEV, TWEETER_AO_CHAN,
			 rangetAo, AO_AREF, tweeterDataMin);
  if (rc<0){
    comedi_perror("comedi_write_data 3");
    exit(1);
  }
  //--- wait 25ms
  usleep(25000);
  }

  //--- set tweeter at 0V
  rc = comedi_data_write(cifh, AO_SUBDEV, TWEETER_AO_CHAN,
			 rangetAo, AO_AREF, tweeterDataZero);
  if (rc<0){
    comedi_perror("comedi_write_data 4");
    exit(1);
  }

#else

  //--- start comedi Analog Output scan
  rc = comedi_command(cifh, &cmdAo);
  fprintf(stderr,"Start AO wave, rc= %d\n", rc);
  if (rc<0) {
    fprintf(stderr,"errno=%d\n",errno);
    comedi_perror("comedi_command");
    exit(1);
  }

  unsigned int aoScanDataCycle(0);
  rc = comedi_write_wave(cifh,
			 aoScanData, sizeof(aoScanData)/sizeof(sampl_t),
			 aoScanDataCycle, AOSCAN_SAMPLES);
  if (rc<0) {
    perror("comedi_write_wave");
    exit(1);
  }
  
  rc = comedi_internal_trigger(cifh, AO_SUBDEV, 0);
  if (rc<0) {
    perror("comedi_internal_trigger");
    exit(1);
  }
  
#ifdef CONTINUOUS_TEST
  for(;;) {
    rc = comedi_write_wave(cifh,
			   aoScanData, sizeof(aoScanData)/sizeof(sampl_t),
			   aoScanDataCycle, AOSCAN_SAMPLES);
  }
#endif

#endif

  //
  // initialize Root
  //
  // code inspired from:
  // - test/hsimple.cxx of Root 3.05.3 package
  //

  // Create a new ROOT binary machine independent file.
  // This file is now becoming the current directory.
  TFile hfile("labFisOndeSup.root", "RECREATE",
	      "Lab Fis: Onde Superficiali");

  //--- allocate Root Tree
  TTree *tree = new TTree("nt", "Lab Fisica: Onde Superficiali");

  //--- Root ntuple definition: see Root/tutorials/cernbuild.C

  //--- declare & alloc floating point vars for photo detector voltages
  Float_t phDet1;
  Float_t phDet2;
  //--- add to "tree" ntuple branches for photo detector voltages
  tree->Branch("phDet1",&phDet1,"phDet1/F");
  tree->Branch("phDet2",&phDet2,"phDet2/F");

  TH1F *hsPhDet1 = new TH1F("hsPhDet1", "Photodetector #1", 100, -5, 5);
  hsPhDet1->SetXTitle("Volts");
  hsPhDet1->SetYTitle("events");

  TH1F *hsPhDet2 = new TH1F("hsPhDet2", "Photodetector #2", 100, -5, 5);
  hsPhDet1->SetXTitle("Volts");
  hsPhDet1->SetYTitle("events");

  TProfile *hsPhDet1TimeProf =
    new TProfile("hsPhDet1TimeProf", "Photodetector #1 vs. time",
		 AISCAN_SAMPLES/20, -0.5/AISCAN_RATE,
		 (static_cast<double>(AISCAN_SAMPLES)-0.5)/AISCAN_RATE,
		 -10, 10);
  hsPhDet1TimeProf->SetXTitle("seconds");
  hsPhDet1TimeProf->SetYTitle("Volts");

  TProfile *hsPhDet2TimeProf =
    new TProfile("hsPhDet2TimeProf", "Photodetector #2 vs. time",
		 AISCAN_SAMPLES/20, -0.5/AISCAN_RATE,
		 (static_cast<double>(AISCAN_SAMPLES)-0.5)/AISCAN_RATE,
		 -10, 10);
  hsPhDet2TimeProf->SetXTitle("seconds");
  hsPhDet2TimeProf->SetYTitle("Volts");

  //--- ask for non-blocking I/O on comedi
  fcntl(comedi_fileno(cifh),F_SETFL,O_NONBLOCK);

  //--- allocate object to setup terminal input
  labFisTermIoHandler tioh(stdin);
  //--- setup terminal for immediate keypress input with no echo
  tioh.SetImmedNoecho();

  //--- init sample number
  scansTot = 0;

  //--- DAQ loop (endless)
  while (1) {
    //--- read available data
    rc = labFisScanRead(cifh, AI_SUBDEV, &aiScanData[0], sizeof(aiScanData),
			static_cast<long int>(1e6), 10);

    //--- end endless loop if end of data reached
    if (rc == 0) {
      if (pendingNewline) fprintf(stderr, "\n");
      // fprintf(stderr, "EOF, ending scan after reading %6d samples\n",
      //     scansTot);
      break;
    }

    //--- read available data
    if (rc > 0) {
      unsigned int i;

      len = rc/sizeof(sampl_t);
      for(i=0; i < len; i++, scansTot++) {

	unsigned int scanNum = scansTot / AISCAN_NUM;
	Float_t sampleTime = static_cast<double>(scanNum)/AISCAN_RATE;

	switch(scansTot % AISCAN_NUM) {
	case 0:
	  phDet1 = comedi_to_phys(aiScanData[i], rangepAi, maxDataAi);
	  //--- fill histogram of read voltages
	  hsPhDet1->Fill(phDet1);
	  //--- fill histogram of read voltages vs. time in seconds
	  hsPhDet1TimeProf->Fill(sampleTime, phDet1);
	  break;
	case 1:
	  phDet2 = comedi_to_phys(aiScanData[i], rangepAi, maxDataAi);
	  //--- fill histogram of read voltages
	  hsPhDet2->Fill(phDet2);
	  //--- fill histogram of read voltages vs. time in seconds
	  hsPhDet2TimeProf->Fill(sampleTime, phDet2);

	  //--- once the 2nd and last element is filled, save the ntuple
	  tree->Fill();
	  break;
	}

	// fprintf(stderr, "%4x ", aiScanData[i]);
	pendingNewline = 1;
	if (i%16==15) {
	  // fprintf(stderr,"\n");
	  pendingNewline = 0;
	}
      }
      continue;
    }

    if (-rc == LFDAQ_STDINRDY) {
      unsigned char ch = fgetc(stdin);
      if (ch == 25) {
	fprintf(stderr,"  --- <control-Y> received, quitting ... --- \n");
	break;
      } else {
	fprintf(stderr,"  --- Press <control-Y> to quit --- \n");
      }
      continue;
    }

    if (-rc == LFDAQ_TIMEOUT) {
      fprintf(stderr,"Timeout waiting for data \n");
      break;
    }

    break;
  }

  //--- Save all objects in this file
  hfile.Write();

  //--- Close the file.
  hfile.Close();

  //--- open saved file for reading
  TFile hfileR("labFisOndeSup.root", "READ");

  Float_t hMax, hMin;

  //--- recover a pointer to the histogram "hsPhDet1TimeProf"
  TH1D *hs1 = (TH1D *)hfileR.Get("hsPhDet1TimeProf");
  hMax = hs1->GetMaximum();
  hMin = hs1->GetMinimum();
  hs1->SetMaximum(hMax+0.1*(hMax-hMin));
  hs1->SetMinimum(hMin-0.1*(hMax-hMin));

  //--- recover a pointer to the histogram "hsPhDet2TimeProf"
  TH1D *hs2 = (TH1D *)hfileR.Get("hsPhDet2TimeProf");
  hMax = hs2->GetMaximum();
  hMin = hs2->GetMinimum();
  hs2->SetMaximum(hMax+0.1*(hMax-hMin));
  hs2->SetMinimum(hMin-0.1*(hMax-hMin));

  //--- init an object to get Root GUI on a canvas
  TApplication theApp("App", &argc, argv);

  //--- init an object to get a graphical canvas
  TCanvas *canvas = new TCanvas("canvas", "Photo detectors vs. time",
				800, 800);
  canvas->Divide(1,2);

  //--- draw the histograms
  canvas->cd(1);
  hs1->Draw("hist C");
  canvas->cd(2);
  hs2->Draw("hist C");

  //--- get the canvas updated
  canvas->Update();

  //--- recover pointer to Root tree (ntuple) "nt"
  // TTree *treeR = (TTree *)hfileR.Get("nt");

  //--- start root GUI: "Exit ROOT" will return to the next statement
  theApp.Run(kTRUE);

  //--- Close the file.
  hfile.Close();

  return 0;
}
